- Plot vector geospatial data (points/polygons)
- Plot raster geospatial data (image)
- Extracting data from raster
- Making maps (base maps, compass rose, overlays, etc)
2020-02-26
We are researchers who want to answer the question "What range of environmental conditions do zebra tend to inhabit?"
General workflow:
Variants of this question/workflow are applicable to other study systems.
sf package for vector data
sp, rdgal for read/write, and rgeos for geometric operationstidyverseraster package for raster operations and more
sf objectssf (simple feature) objects are extended data.frames or tibbles
tibblesfc (simple feature columns) with useful components like
bboxCRS (coordinate reference system) attributessf and tidyversesf functions begin with st_summary, plotsf methods for filter, arrange, distinct, group_by, ungroup, mutatemethods(class = "sf") ## [1] $<- [ ## [3] [[<- aggregate ## [5] anti_join arrange ## [7] as.data.frame cbind ## [9] coerce dbDataType ## [11] dbWriteTable df_spatial ## [13] distinct extent ## [15] extract filter ## [17] full_join gather ## [19] group_by group_map ## [21] group_split identify ## [23] initialize inner_join ## [25] left_join mask ## [27] merge mutate ## [29] nest plot ## [31] print raster ## [33] rasterize rbind ## [35] rename right_join ## [37] sample_frac sample_n ## [39] select semi_join ## [41] separate separate_rows ## [43] show slice ## [45] slotsFromS3 spread ## [47] st_agr st_agr<- ## [49] st_area st_as_sf ## [51] st_bbox st_boundary ## [53] st_buffer st_cast ## [55] st_centroid st_collection_extract ## [57] st_convex_hull st_coordinates ## [59] st_crop st_crs ## [61] st_crs<- st_difference ## [63] st_force_polygon_cw st_geometry ## [65] st_geometry<- st_interpolate_aw ## [67] st_intersection st_intersects ## [69] st_is st_is_polygon_cw ## [71] st_join st_line_merge ## [73] st_linesubstring st_make_valid ## [75] st_minimum_bounding_circle st_nearest_points ## [77] st_node st_normalize ## [79] st_point_on_surface st_polygonize ## [81] st_precision st_segmentize ## [83] st_set_precision st_simplify ## [85] st_snap st_snap_to_grid ## [87] st_split st_subdivide ## [89] st_sym_difference st_transform ## [91] st_transform_proj st_triangulate ## [93] st_union st_voronoi ## [95] st_wrap_dateline st_write ## [97] st_zm summarise ## [99] transmute ungroup ## [101] unite unnest ## see '?methods' for accessing help and source code
roads <- st_read(here("data", "enproads.shp"))
## Reading layer `enproads' from data source `C:\Users\mbuon\Desktop\QERM597-making-maps\data\enproads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 565 features and 3 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 14.35337 ymin: -19.48418 xmax: 17.15124 ymax: -18.5008
## epsg (SRID): NA
## proj4string: NA
class(roads)
## [1] "sf" "data.frame"
head(roads)
## Simple feature collection with 6 features and 3 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 14.46841 ymin: -19.41752 xmax: 14.58964 ymax: -19.33077
## epsg (SRID): NA
## proj4string: NA
## TYPE LENGTH KM geometry
## 1 Track 0.7 1 LINESTRING (14.50502 -19.35...
## 2 Track 3.6 4 LINESTRING (14.49077 -19.36...
## 3 Graded 2.3 2 LINESTRING (14.4864 -19.417...
## 4 Track 7.1 7 LINESTRING (14.58498 -19.34...
## 5 Track 1.2 1 LINESTRING (14.58498 -19.34...
## 6 Graded 0.1 0 LINESTRING (14.52596 -19.33...
# Highlight parts
# 1. list column (aka sfc)
# 2. feature geometry (sfg)
# 3. feature (row)
# Default methods for objects plot(roads)
plot(filter(roads[,3], KM > 20))
rnaturalearthSpatial* objects using st_as_sfworld <- ne_countries(scale = "medium", returnclass = "sf") st_geometry(world) %>% plot()
Default data is Lat/Long
Namibia <- ne_countries(country = "namibia") Namibia_sf <- st_as_sf(Namibia) plot(Namibia_sf)
* Exercise: try for another country
#Uncomment here if using rnatural earth for the first time:
# devtools::install_github("ropenscilabs/rnaturalearthhires", force = TRUE)
# install.packages("rnaturalearthhires", repos = "http://packages.ropensci.org")
namibia_political <- ne_states(country = "namibia")
namibia_political_sf <- st_as_sf(namibia_political)
plot(namibia_political_sf)
st_crs("+proj=longlat +datum=WGS84") # Proj.4 string"
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(3857) # ESPG = 3857
## Coordinate Reference System:
## EPSG: 3857
## proj4string: "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
st_crs(3857)$units # check units
## [1] "m"
st_crs(4326)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(4326)$units
## NULL
st_crs(NA) # unknown (assumed planar/Cartesian)
## Coordinate Reference System: NA
st_transform transforms or converts coordinates to new reference system(a1 <- st_area(Namibia_sf)) ## 824886560992 [m^2] (a2 <- st_area(st_transform(Namibia_sf, 32733))) # Namibia plane m (UTM) ## 826175760745 [m^2] (a3 <- st_area(st_transform(Namibia_sf, 4293))) # Another reference from espg.io ## 824892086801 [m^2] (a4 <- st_area(st_transform(Namibia_sf, 4326))) ## 824886560992 [m^2] # Namibia projection = 4326
*If you are working with U.S. data you can convert to U.S. feet
Change the units!
units::set_units(a1, km^2) ## 824886.6 [km^2] units::set_units(a1, ft^2) ## 8.879005e+12 [ft^2]
zebra_data <- read.csv(here("data", "zebra.csv"))
head(zebra_data)
## X.1 X Name Date Longitude Latitude Speed Direction Temp
## 1 50 78 AG059 4/20/2009 15:40 15.89605 -19.12138 0 0 0
## 2 34 60 AG059 4/20/2009 13:01 15.91854 -19.17761 0 0 0
## 3 33 59 AG059 4/20/2009 12:02 15.91857 -19.17766 0 0 0
## 4 62 90 AG059 4/20/2009 17:00 15.80543 -19.08002 0 0 0
## 5 109 159 AG059 4/21/2009 1:40 15.91887 -19.17758 0 270 0
## 6 107 155 AG059 4/21/2009 1:00 15.91882 -19.17760 0 270 0
## Altitude PDOP ID DatePOS diff day
## 1 1123 2 AG059 4/20/2009 15:40 19 4/20/2009
## 2 1113 3 AG059 4/20/2009 13:01 59 4/20/2009
## 3 1130 3 AG059 4/20/2009 12:02 721 4/20/2009
## 4 1139 2 AG059 4/20/2009 17:00 17 4/20/2009
## 5 1133 2 AG059 4/21/2009 1:40 40 4/21/2009
## 6 1126 4 AG059 4/21/2009 1:00 20 4/21/2009
zebra_sf <- read.csv(here("data", "Zebra.csv")) %>%
dplyr::select(ID = Name, 4:6) %>%
mutate(timestamp = as.POSIXct(lubridate::mdy_hm(Date))) %>%
st_as_sf(., coords = 3:4, crs = "+init=epsg:4326")
st_transform(zebra_sf, 32733) #convert to UTM
## Simple feature collection with 11490 features and 3 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 546746.9 ymin: 7852220 xmax: 691330.5 ymax: 7912185
## epsg (SRID): 32733
## proj4string: +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs
## First 10 features:
## ID Date timestamp geometry
## 1 AG059 4/20/2009 15:40 2009-04-20 15:40:00 POINT (594243.6 7885501)
## 2 AG059 4/20/2009 13:01 2009-04-20 13:01:00 POINT (596576 7879266)
## 3 AG059 4/20/2009 12:02 2009-04-20 12:02:00 POINT (596579.5 7879260)
## 4 AG059 4/20/2009 17:00 2009-04-20 17:00:00 POINT (584733 7890124)
## 5 AG059 4/21/2009 1:40 2009-04-21 01:40:00 POINT (596611.3 7879269)
## 6 AG059 4/21/2009 1:00 2009-04-21 01:00:00 POINT (596605.8 7879266)
## 7 AG059 4/21/2009 5:01 2009-04-21 05:01:00 POINT (596604.9 7879271)
## 8 AG059 4/21/2009 4:41 2009-04-21 04:41:00 POINT (596610.9 7879268)
## 9 AG059 4/22/2009 4:00 2009-04-22 04:00:00 POINT (596603.6 7879269)
## 10 AG059 4/22/2009 2:01 2009-04-22 02:01:00 POINT (596609.4 7879272)
countries <- ne_countries(scale=110) p <- ggplot() + geom_polygon(data=countries, aes(x=long, y=lat, group=group), color="white", lwd = .25) ggplotly(p)
st_writeZebra.csv and plot the geometry.roads <- st_read(here("data", "enproads.shp"), crs = "+init=epsg:4326") %>% #4326= coordinate system based on Earth's center of mass
st_transform(32733) #32733 = spatial reference for Nambia
## Reading layer `enproads' from data source `C:\Users\mbuon\Desktop\QERM597-making-maps\data\enproads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 565 features and 3 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 14.35337 ymin: -19.48418 xmax: 17.15124 ymax: -18.5008
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
st_geometry(roads) %>% plot
zebra_sf <- read.csv(here("data", "Zebra.csv")) %>%
dplyr::select(ID = Name, 4:6) %>%
mutate(timestamp = as.POSIXct(lubridate::mdy_hm(Date))) %>%
st_as_sf(., coords = 3:4, crs = "+init=epsg:4326") %>% #longlat, converting foreign object to SF object
st_transform(32733) #convert to UTM
ggplot() +
geom_sf(data=roads) +
geom_sf(data=zebra_sf, aes(color=ID))
How close or far from different road types do zebra move?
unique(roads$TYPE)
## [1] Track Graded Gravel Tar
## Levels: Graded Gravel Tar Track
#Filter roads based off type:
large_roads <- filter(roads, TYPE %in% c("Tar", "Gravel"))
small_roads <- filter(roads, TYPE %in% c("Graded", "Track"))
ggplot() +
geom_sf(data=large_roads, size=1.5) +
geom_sf(data=small_roads, size=0.6) +
geom_sf(data=zebra_sf, aes(color=ID))
# Find the minimum distance (in meters) of each point to a large road large_dist<- st_distance(y=zebra_sf, x=large_roads) # a units matrix dim =[nrow(x), nrow(y)]; takes 10-20 seconds zebra_sf$large_road_dist <- apply(large_dist, 2, min) # Find the minimum distance (in meters) of each point to a small road small_dist<- st_distance(y=zebra_sf, x=small_roads) # a units matrix dim =[nrow(x), nrow(y)]; takes 10-20 seconds zebra_sf$small_road_dist <- apply(small_dist, 2, min) head(data.frame(zebra_sf)) ## ID Date timestamp geometry ## 1 AG059 4/20/2009 15:40 2009-04-20 15:40:00 POINT (594243.6 7885501) ## 2 AG059 4/20/2009 13:01 2009-04-20 13:01:00 POINT (596576 7879266) ## 3 AG059 4/20/2009 12:02 2009-04-20 12:02:00 POINT (596579.5 7879260) ## 4 AG059 4/20/2009 17:00 2009-04-20 17:00:00 POINT (584733 7890124) ## 5 AG059 4/21/2009 1:40 2009-04-21 01:40:00 POINT (596611.3 7879269) ## 6 AG059 4/21/2009 1:00 2009-04-21 01:00:00 POINT (596605.8 7879266) ## large_road_dist small_road_dist ## 1 7.845860 3479.8987 ## 2 163.187201 241.8660 ## 3 156.746071 245.8446 ## 4 9.115608 1605.2154 ## 5 151.769751 227.9759 ## 6 151.358102 231.6984
ggplot(zebra_sf) +
geom_histogram(aes(large_road_dist, fill="large roads"), alpha=0.5) +
geom_histogram(aes(small_road_dist, fill="small roads"), alpha=0.5) +
scale_fill_manual(labels=c("large roads", "small roads"), values=c("blue", "orange"))
Assume you have an existing 'matrix' object in your workspace that you want to cast as a 'raster' object.
wspd.raster<-raster(wspd,xmn=x.min,ymn=y.min,xmx=x.max,ymx=y.max) class(wspd) ## [1] "matrix" class(wspd.raster) ## [1] "RasterLayer" ## attr(,"package") ## [1] "raster"
plot(wspd.raster)
You should have a file labeled "ndvi_mean_utm.tif" in your working directory. Use ?raster and upload the raster into your workspace. Plot a simple raster map.
Solution:
NDVI_mean <- raster(here("data", "ndvi_mean_utm.tif"))
NDVI_mean
## class : RasterLayer
## dimensions : 474, 1274, 603876 (nrow, ncol, ncell)
## resolution : 231.6564, 231.6564 (x, y)
## extent : 431873.8, 727004, 7844489, 7954294 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : C:/Users/mbuon/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif
## names : ndvi_mean_utm
## values : -558.117, 5303.021 (min, max)
Solution:
plot(NDVI_mean)
ggplot(data.frame(rasterToPoints(wspd.raster))) + geom_raster(aes(x=x, y=y,fill=layer))
Take the raster object and plot it in ggplot. Can you overlay the zebra location data from earlier? (EC) Can you write equivalent code using the pipe (i.e. %>%) functionality?
Solution: (no piping)
ggplot(data.frame(rasterToPoints(NDVI_mean))) + geom_raster(aes(x=x, y=y,fill=ndvi_mean_utm)) + geom_sf(data=zebra_sf, size=0.7,col="black")
Solution: (piping)
NDVI_mean %>% rasterToPoints() %>% data.frame() %>% ggplot()
Solution: (no piping)
ggplot(data.frame(rasterToPoints(NDVI_mean)))
pop_sf <- read_csv(here("data", "huc_whp_burnable_area_population.csv")) %>%
st_as_sf(., coords = 6:5, crs = "+init=epsg:4326")
## Parsed with column specification:
## cols(
## huc_name = col_character(),
## whp = col_double(),
## acres_burnable = col_double(),
## population_in_huc = col_double(),
## Latitude = col_double(),
## Longitude = col_double()
## )
pop_sf$wspd <- raster::extract(wspd.raster, as(pop_sf, "Spatial"))
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 2.130 2.927 3.361 3.442 3.974 5.414
Modify this code to extract NDVI information at the Zebra locations. Calculate summary statistics.
Solution:
zebra_sf$ndvi<-raster::extract(NDVI_mean, as(zebra_sf, "Spatial")) summary(zebra_sf$ndvi) ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 626.5 2527.6 2682.2 2775.8 3015.1 4212.1
slice<-raster(ncol=20,nrow=20)
slice[]<-rnorm(n=ncell(slice))
sim_stack <- stack(slice,slice^2,1/slice)
dates<-c(1:3)
#assign time value (in this case year) to each layer
sim_stack <- setZ(sim_stack, dates)
#time match and extract
stack_dates <- getZ(sim_stack) #all raster layer dates
zebra_sf$NDVI <- NA
for (i in 1:50) { #this would take a long time to run through nrow(zebra_sf) locations, so demonstrating for first 50 locations
print(paste0("extracting zebra location #",i))
zeb_date <- as.Date(zebra_sf$timestamp[i])
stack_idx <- which(stack_dates %in% seq(zeb_date - 15, zeb_date, by='day')) ##which raster layer date matches zebra date
zebra_sf$NDVI[i] <-raster::extract(NDVI_stack[[stack_idx]], #as(zebra_sf$geometry[i], "Spatial")) #extract
}
head(data.frame(zebra_sf))
Using netCDF files (ncdf4 package in R) are my preferred approach.
Focus on ggmap
ggmap GitHub page for more information about API keysggmapggmapTwo steps:
This is my personal API key, which you can use for the purposes of this class!
register_google(key = "AIzaSyBRkw_wzDKuWZDikdD86Kmp1Sa6PzuCKFc")
To add a Stamen basemap, first define the bounding box for the basemap you would like to download.
zebra_bbox <- c(14, -20, 17.5, -18)
names(zebra_bbox) <- c("left", "bottom", "right", "top")
terr_basemap <- get_stamenmap(bbox=zebra_bbox, zoom=9, maptype="terrain") ggmap(terr_basemap) + geom_sf(data=roads, inherit.aes = FALSE) + geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) + coord_sf(crs = st_crs(4326))
wat_basemap <- get_stamenmap(bbox=zebra_bbox, zoom=9, maptype="watercolor") ggmap(wat_basemap) + geom_sf(data=roads, inherit.aes = FALSE) + geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) + coord_sf(crs = st_crs(4326))
To add a Google basemap, first define the center coordinates for the basemap you would like to download.
zebra_center <- c(15.8, -19)
names(zebra_center) <- c("lon", "lat")
sat_basemap <- get_googlemap(center=zebra_center, zoom=8, size=c(640, 640), scale=2,
maptype="satellite")
ggmap(sat_basemap) +
geom_sf(data=roads, inherit.aes = FALSE) +
geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326))
Notice that get_googlemap() returns a square basemap tile, which then sets the plotting limits of your map.
To manually set different plotting limits, pass in an empty "base layer" and set x and y limits using ggplot().
ggmap(sat_basemap, maprange=FALSE, base_layer=ggplot(aes(x=1, y=1), data=NULL)) +
xlim(14.3, 17.4) + ylim(-20, -18) + xlab("lon") + ylab("lat") +
geom_sf(data=roads, inherit.aes = FALSE) +
geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326))
RgoogleMaps
ggplot2 (boo)tif file as a raster brick, then display red/green/blue valuesgeocode function to find the lat/lon coordinates for a known location
geocode("Space Needle")my_basemap <- get_googlemap(center="Space Needle", zoom=14, size=c(640, 640),
scale=2, maptype="satellite")
ggmap(my_basemap)
q0 <- ggplot() +
geom_sf(data = large_roads, size = 1.5, color = "grey", inherit.aes = FALSE) +
geom_sf(data = small_roads, size = 0.6, color = "grey", inherit.aes = FALSE) +
geom_sf(data = zebra_sf, aes(color = ID), inherit.aes = FALSE,
alpha = 0.1, # fix over plotting
show.legend = "point") # make legend keys points instead of rects.
q1 <- q0 +
coord_sf(xlim = c(14, 17.5), # make space for legend
ylim = c(-20.5, -18.25),
expand = FALSE,
crs = st_crs(4326)) +
scale_color_manual(
values = viridis::viridis(10), # pick new color scale
name = "Zebra ID", # change title of legend
guide = guide_legend(
direction = "horizontal", # horizontal legend
title.position = "top", # title at top
title.hjust = 0.5, # center title
byrow = TRUE, # organize legend keys
nrow = 3, # organize legend keys
override.aes = list(alpha = 1) # make the legend keys 100% opaque
)) +
theme_bw() +
theme(legend.position = c(0.65, 0.2), # position legend in lower R corner
legend.key = element_rect(fill = NA), # remove background of legend keys
panel.grid.major = element_line(color = "white")) + # remove graticules
labs(x = "Longitude", y = "Lattitude")
q2 <- q1 +
annotation_scale(location = "bl") + # from ggspatial package
annotation_north_arrow(location = "bl", which_north = "true",
pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
style = north_arrow_fancy_orienteering)
ne_namibia <- ne_countries(scale = "medium", country = "namibia", returnclass = "sf")
n1 <- ggplot(ne_namibia) + # need to put data in ggplot call for geom_rect
geom_sf() +
geom_sf(data = large_roads, inherit.aes = FALSE) +
geom_sf(data = small_roads, inherit.aes = FALSE) +
geom_rect(xmin = 14, xmax = 18, ymin = -20, ymax = -18, # extent indicator
fill = NA, colour = "black", size = 1.5) +
theme_bw() +
coord_sf(xlim = c(11, 47),
ylim = c(-34, -16.9), # make additional space for inset map
expand = FALSE,
datum = NA) # remove graticules and coords
# Calculate ratio of inset plot
# use st_bbox(zebra_sf)
inset_ratio <- (-18.5 - (-20.5))/(17.5 - 14)
myfig <- cowplot::ggdraw(n1) +
cowplot::draw_plot(q2, width = 1.2, height = 1.2*inset_ratio,
x = 0.01, y = 0.1, hjust = 0, vjust = 0)